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Noise and Correlations in a Spatial Population Model with Cyclic Competition 
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Noise and spatial degrees of freedom characterize most ecosystems. Some aspects of their influ- 
ence on the coevolution of populations with cyclic interspecies competition have been demonstrated 
in recent experiments [e.g. B. Kerr et al., Nature 418, 171 (2002)]. To reach a better theoretical 
understanding of these phenomena, we consider a paradigmatic spatial model where three species 
exhibit cyclic dominance. Using an individual-based description, as well as stochastic partial dif- 
ferential and deterministic reaction-diffusion equations, we account for stochastic fluctuations and 
spatial diffusion at different levels, and show how fascinating patterns of entangled spirals emerge. 
We rationalize our analysis by computing the spatio-temporal correlation functions and provide 
analytical expressions for the front velocity and the wavelength of the propagating spiral waves. 

PACS numbers: 87.23.Cc,02.50.Ey,05.10.Gg,87.18.Hf 



Understanding the combined influence of spatial de- 
grees of freedom and noise on biodiversity is an impor- 
tant issue in theoretical biology and ecology. This im- 
plies to face the challenging problem of studying com- 
plex noncquilibrium structures, which form in the course 
of nonlinear evolution 0, 0, H, 0, IS E|- More gener- 
ally, self-organized nonequilibrium patterns and travel- 
ing waves are ubiquitous in nature and appear, for in- 
stance, in chemical reactions, biological systems, as well 
as in epidemic outbreaks [7|. Among the most stud- 
ied types of patterns are spiral waves, which are rel- 
evant to autocatalytic chemical reactions, aggregating 
slime-mold cells and cardiac muscle tissue [8|. In all 
these nonequilibrium and norzimearprocesses, as well as 
in population dynamics models 0, [j, HI , pattern forma- 
tion is driven by diffusion which, together with internal 
noise, act as mechanisms allowing for stabilization and 
coevolution of the reactants. In this work, we consider 
a paradigmatic spatially-extended 3 species population 
system with cyclic competition, which can be regarded 
as a simple food-chain model [9(. In fact, such a system 
is inspired by recent experiments on the coevolution of 3 
species of bacteria in cyclic competition Q . Using meth- 
ods of statistical physics, wc study the influence of spatial 
degrees of freedom and internal noise on the coevolution 
of the species and on the emerging spiral patterns. In par- 
ticular, we compute the correlation functions and provide 
analytical expressions for the spreading speed and wave- 
length of the propagating fronts. To underpin the role 
of internal noise, the results of the stochastic description 
are compared with those of the deterministic equations. 

In this Letter, we investigate a stochastic spatial vari- 
ant of the rock-paper- scissors game [j| (also referred to 
as cyclic Lotka-Volterra model). These kinds of systems 
have been studied both from a game-theoretic perspec- 



tive, see e.g. |1(X 1 1 11 ) and references therein, and within 
the framework of chemical reactions [!, HI] , revealing rich 
spatio-temporal behaviors (e.g. emergence of rotating 



spirals). While our methods have a broad range of ap- 
plicability, they are illustrated for a prototypical model 
introduced by May and Leonard [13j where 3 species, 
A, B and C undergo a cyclic competition (codominancc 
with rate a) and reproduction (with rate /i), according 
to the reactions 
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Hence, an individual of species A will consume one of 
species B {AB — ► A(d) with rate a and will reproduce 
with rate p if an empty spot, denoted 0, is available 
(A<Z) —* AA, i.e. there is a finite carrying capacity). In 
addition, to mimic the possibility of migration, it is re- 
alistic to endow the individuals with a form of mobility. 
For the sake of simplicity, we consider a simple exchange 
process, with rate e, among any nearest-neighbor pairs 
of agents: XY YX, where X,Y £ {A,B,C,<Z>}. If 
one ignores the spatial structure and assumes the system 
to be well- mixed (with an infinite number of individuals), 
the population's mobility plays no role and the dynam- 
ics is aptly described by the deterministic rate equations 
(RE) for the densities a, b, c of species A, B and C, re- 
spectively. Introducing s = (a,b,c), the RE read: 



d t Si = Si[fi(l - p) - o-s i+2 ], i G {1, 2, 3} 



(2) 



where the index i is taken modulo 3 and p = a + b + c 
is the total density. As shown by May and Leonard [HI 
(see also [Hj]), these equations possess 4 absorbing fixed 
points, corresponding to a system filled with only one 
species and to an empty system. In addition, there is 
a reactive fixed point s* = a + 3 (1, 1, 1), corresponding 
to a total density p* = a ^ ■ A linear stability analysis 
shows that s* is unstable. The absorbing steady states 
(1, 0, 0), (0, 1, 0) and (0, 0, 1) are heteroclinic points. The 
existence of a Lyapunov function C = abc/ p 3 allows to 
prove that, within the realm of the above RE, the phase 
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portrait is characterized by flows spiraling outward from 
s*, with frequency u>o — V3pcr/[2(3p+2a)] in its vicinity. 
Approaching the boundaries of the phase portrait, the 
trajectories form (heteroclinic) cycles indefinitely close 
to the edges (without ever reaching them), with densi- 
ties approaching in turn the value one. Despite its math- 
ematical elegance, this behavior has been recognized to 
be unrealistic 0, EH- I n iac tj f° r finite populations, 
fluctuations arise and always cause the extinction of two 
species in finite time (see e.g. Ref [HI). 

In this work, considering the spatial version of the 
above model in the presence of internal noise, wc show 
that a robust (and, arguably, more realistic) scenario 
for the evolution arises. The reaction schemes ([T]) 
and the exchange events are considered to occur on a 
d— dimensional regular lattice of N sites, labeled r = 
(ri, . . . ,rd). Each lattice site has z neighbors at a dis- 
tance 5r (e.g. z = 2d and N = L d for hypercubic lat- 
tices of linear size L) and is either empty or occupied 
by at most one individual. On the lattice, the binary 
reactions (JlJ and exchanges only occur among pairs of 
nearest-neighbors. In the situation of large system sizes, 
the continuum limit reveals that for the exchange pro- 
cess to be an efficient driving mechanism, the rate e has 
to scale as £ a N v , with v = 2/d and N — > oo. In fact, 
if < v < 2/d the system is dominated by the local 
reactions ([1} among neighboring individuals; while effec- 
tive diffusion renders locality irrelevant when v > 2/d. 
Only when v = 2/d, there is an effective competition be- 
tween the stirring process and the local reactions fl}. It 
is therefore useful to introduce the effective diffusion con- 
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number of individuals involved in the reactions, internal 
fluctuations arise in the system. The latter originate from 
(i) the interspecies reactions (JT]) and (ii) the exchange 
processes. In the continuum limit, where Sr = N^ x / d 
with TV, e — > oo (and finite D), there is a separation of 
time scales and the pair exchanges occur much faster 
than the reactions (e oc N 2 / d ). Actually, the fluctua- 
tions associated with ([T|) and the agents' mobility scale 
respectively as TV -1 / 2 and TV -1 , with the former domi- 
nating over the latter and being the only relevant con- 
tribution. This result is revealed by a system size, also 
called Kramers- Moyal (see e.g. Ref. [16(, Chap. 8), ex- 
pansion (SZE) of the master equation underlying the ex- 
change processes and the reactions (JXJ) [T3] • Furthermore, 
the SZE yields a proper Fokkcr-Planck equation, which is 
equivalent to a set of (Ito) stochastic partial differential 
equations (SPDE) with white noise. The derivation, ob- 
tained in the continuum limit from the master equation, 
is outlined in the supplementary EPAPS document 
and will be detailed elsewhere [l8j. Here, we quote the 
expression of the SPDE: 




FIG. 1: (Color online) Snapshots of reactive steady states for 
rates D = 3x10~ 6 ,/i = <j = 1. Each color (level of grey) 
indicates one species (black dots correspond to vacancies). In 
(a) and (b) results are from lattice simulations for L = 200 (a) 
and L — 1000 (b), i.e. different e. Spiral structures emerge 
for sufficiently large exchange rate (b). Numerical solution 
of the SPDE and DRDE are shown in (c), resp. (d); see 
text. In (a), (b) and (c), initially s(r,0) = s* . 



where V 2 is the Laplacian operator; (£,i(r,t)) = 0, 
(&(r,%(r',t')) = S itj S(r - r')6(t - f) and 



Ai(s) = Si\p(l - p) -as i+2 ], (4) 
Cij(s) = S ij y/N- 1 s i [fj,(l-p) + cTSi + 2\. (5) 

Again, the indices are taken modulo 3 and now Sj = 
Si(r,t). As explained in 1(| 3|, these SPDE have to be 
interpreted in the sense of Ito calculus. While Eqs . ([3]) 



d tSi =DV 2 s i +A l (s) 



E 



Cij(*)& *€ {1,2,3} (3) 



and our approach are valid in any dimension 17|, for 
specificity, we now analyze the spatio-temporal proper- 
ties of the system in two dimensions with periodic bound- 
ary conditions. On the one hand we have solved numer- 
ically the SPD E dp using the open software from the 
XMDS project [19(. On the other hand, we have car- 
ried out individual-based simulations of the reactions ((T|) 
for mobile (exchange process) particles on lattices of size 
Lx L, with L = 30—1000. This allows to carefully study 
the convergence towards the continuum limit, where the 
description in terms of ([3]) is expected to be accurate. 

As other spatially-extended dynamical systems 0, 
11 . l^ |, the model under consideration displays fascinat- 
ing nonequilibrium patterns emerging in the course of the 
evolution. In Fig. Q] (a) and (b) , we report typical long- 
time snapshots of the system for low (a) and high (b) ex- 
change rates (but keeping D fixed) , as obtained from lat- 
tice simulations. In both cases we notice that intriguing 
patterns form. For slow exchange rate, the system dis- 
plays non-geometrical patches, similarly to what happens 
in systems with self-organized criticality [2(|. When the 
exchange rate is raised, the patterns display spiral struc- 
tures. In fact, starting from a spatially homogeneous 
initial condition, s(r,0) = s*, the system is randomly 
perturbed by the internal noise and the resulting spa- 
tial inhomogeneities grow and form wavefronts moving 
through the system. The emergence of spiral patterns is 
a feature shared by other excitable systems (see e.g. 
and corresponds to the ability of the system to sustain 



3 



0.15 



0.05 




0.1 

(distance) r 



FIG. 2: (Color online) Spatial correlation functions in 2D, 
obtained from lattice simulations (red, circles; e = 6, L — 
1000), from the solution of the SPDE (0 (dark blue, squares) 
and of the DRDE (green, triangles), see text. The reaction 
rates are fj, = a = 1 and D — 3 x 10~ 6 . Inset: The correlation 
length icon, for D — 5 x 10~ 5 ,/i = a = 1, as function of e 
(i.e. for different lattice sizes) compared to the prediction of 
the SPDE (black line). The latter is in excellent agreement 
with lattice simulations already for e > 5 (i.e. L > 225). 



the propagation of oscillating waves. For sufficiently large 
e, one observes a striking resemblance between the size 
and structure of the patterns obtained from the lattice 
simulations [Fig. CD (b)] and those from the SPDE © 
[Fig. Q] (c)]. To further compare the predictions of the 
SPDE with the lattice simulations, and to gain addi- 
tional information on the structure of the the emerging 
patterns, we have computed the correlation functions, 
g SiS .(r-r',t) = (a i (r,t)s i (r',t))-(s i (r,t))(s J -(r',t)) in 
two dimensions. In Fig. [2] (red and blue curves), we re- 
port the results for g aa (r, t) in the steady state and notice 
an excellent agreement between the results of the lattice 
simulations and the predictions of the SPDE ([3]). The 
inset of Fig. [2j displays the correlation length t co „ (2f| 
as function of e (D is kept fixed, L varies) obtained in 
the lattice simulations, which is found to coincide with 
the prediction of the SPDE already for e > 5. We have 
also computed the autocorrelation function g SiSj (0, t) and 
found, both in the lattice simulations and from the solu- 
tions of the SPDE, an oscillating behavior with a similar 
characteristic frequency, markedly different from [3 ■ 
This confirms that, even for finite exchange rates, the 
solution of the SPDE ([3]) provides an excellent approxi- 
mation of the lattice simulations of the system. This is 
rather surprising since Eqs. ([3|) have been derived in the 
continuum limit, where A" and e — * oo. A comparable 
influence of finite exchange rate in a predator-prey sys- 
tem has been reported recently 2l| . According to the 
SPDE j3]), 4orr scales as D 1 / 2 , so that by raising the dif- 
fusion one increases the size of the spirals. As we have 
shown in Ref. jf|, this happens up to a critical value D c 
(e.g. D c w 4.5 ± 0.5 x KT 4 for fi = a = 1): above 
that threshold, the spiral structures outgrow the system 



size and only one species survives, corresponding to an 
absorbing steady state predicted by Eqs. J2]). 

As the properties of the lattice simulations are well 
captured by the SPDE where the strength of the 
noise scales as A^ -1 / 2 , with A" — > oo, it is natural to 
investigate the actual influence of this internal noise on 
the steady state of the system. To address this issue, 
we have solved numerically (in 2D, with periodic bound- 
ary conditions) the deterministic reaction-diffusion equa- 
tion (DRDE) obtained from (J3j) by dropping the noise 
terms, i.e. d t Si = -DV 2 s; + Ai(s). Of course, to ob- 
tain a nontrivial steady state for the DRDE one has to 
assume spatially inhomogeneous initial conditions. In 
Fig. Q] (d), we have reported a snapshot of the long- 
time behavior predicted by the DRDE starting from 
s(r,0) = s* + (jig cos27rrir 2 , 0, 0). In this case, the 
dynamics evolves towards a reactive steady state which 
also exhibits spiral waves. However, the latter do not 
form entangled structures, but ordered geometrical pat- 
terns. As an example, only four spirals cover the system 
in Fig. Q] (d) [while noise leads to 106 entangled spirals 
in Fig. Q] (c)]. The correlation functions associated with 
the DRDE therefore exhibit only weakly damped spa- 
tial oscillations (see Fig. [51 green triangles). By ana- 
lyzing typical snapshots like those of Fig. Q] (d), we have 
noted that in the deterministic and stochastic [i.e. lattice 
simulations with "large" e and solutions of Eqs. {5J] de- 
scriptions, the spiral waves share the same propagation 
velocity, frequency and wavelength. However, a major 
difference between these descriptions lies on the crucial 
dependence of the DRDE on initial conditions, which de- 
termine the overall number of spirals and their size. On 
the contrary, because the internal noise acts a random 
source of spatial inhomogencities, the lattice stochastic 
system and the SPDE display robust features. In par- 
ticular, we have found noise to induce a universal spiral 
density of about 0.5 per square wavelength. 

Analytical expressions for the spreading velocity and 
the wavelength of the propagating fronts of the DRDE 
can be obtained by considering the dynamics on the in- 
variant manifold of the RE [22j, given by M. : {yc = 



q-(3/jH 



4/i(3/i+2er) (MA 
V3 -V3^ 
-1 2 -1 
lily 

der, the DRDE can be recast in the form of a forced 
complex Ginzburg-Landau equation (CGLE) [23|, [24| . 
By performing the nonlinear transformation za = 

3/J+o- r. /o„,2 i in„, .„._ . /o„.2 1 anc [ 
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^VSyAyB — 5yfj] , upon ignoring nonlincar- 
ities like (Vza,b) 2 , one is left with the following CGLE 
in the variable z = za + izb II 81 ] : 



d t z = DV 2 z+( Cl -iujo)z-C2(l + ic3)\z\ 2 z, (6) 

With C l = 2(37^)> C 2 = " (3 5M ) £+2a) 1 ' T) ' and C 3 = 
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FIG. 3: (Color online) Plot of A, where A = A-/D is the 
spirals' wavelength of the propagating spiral waves. Analyt- 
ical results (red curve, rescaled by a factor 1.6; see text) are 
compared with the solution of the SPDE (|3]) (black circles). 



on the diffusivity: above a certain threshold the system 
is covered by one species 0] . In the absence of noise, the 
equations still predict the formation of spiral waves, but 
their spatial arrangement depends on the initial condi- 
tions. 
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^ui^+iia) ■ The general theory of front propagation 23, 
24J predicts that Eq. ^ always admits traveling waves as 
stable solutions (i.e. no Bcnjamin-Feir or Eckhaus insta- 
bilities occur). We have determined such periodic solu- 
tions by computing, from the dispersion relation of (|6j). 
the spreading velocity v and the spirals' wavelength A 
(details will be given in (lij|): 



2VciT>, A = 2ttc 3 



D \ 1 — a/1 



(7) 



In the stochastic version of the model, the wavelength 
and velocity of the wavefronts have been found to agree 
with those of the deterministic treatment. Hence, the ex- 
pressions (J7J) also apply (for large e, with D < D c ) to the 
results of lattice simulations (rescaled by a factor L) and 
to the solution of the SPDE j3]). For instance, on a square 
grid with p, = a = 1, lattice simulations and Eqs. ([3|) 
yield v w 0.63Z) 1 ' 2 L, in good agreement with the predic- 
tion of (J?]): v = (D/2) 1 / 2 L. For the spirals' wavelength, 
numerical results (lattice simulations and SPDE) yield 
A cx D 1 / 2 as predicted by ([Jj. In Fig. [3l the analytical 
prediction ([Jj for A is compared with the values obtained 
from the SPDE ([3]), yielding a remarkable agreement for 
the functional dependence on the parameter fj,. Yet, as 
Eq. ([6]) does not account for all nonlincaritics, the ana- 
lytical and numerical values differ by a prefactor w 1.6 
(considered in Fig. [3|) [18| . It can still be noted that (El) 
and the predictions (UJ are valid in all dimensions 17llla|. 

Motivated by recent experiments 0], we have con- 
sidered a spatially-extended model with three species 
in cyclic competition, and focused on the spatial and 
stochastic effects. The local character of the reactions 
and internal noise allow mobile populations to coexist 
and lead to pattern formation. We have shown that 
already for finite mobility the lattice model can be de- 
scribed by SPDE. With the latter and lattice simula- 
tions, we have studied how entanglement of spirals form 
and obtained expressions for their spreading velocity and 
wavelength. The size of the patterns crucially depends 
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Supplementary Material (EPAPS Document) 

http://netserver.aip. org/cgi-bin/epaps?ID=E-PRLTAO-99-046747 

FROM THE MASTER EQUATION TO STOCHASTIC PARTIAL DIFFERENTIAL EQUATIONS 

In this Supplementary Notes, we want to explicitly outline how, starting from the master equation associated with 
the stochastic May-Leonard model [defined by the reactions (1)], the set of stochastic partial differential equations 
(3) can be obtained via system size expansion. 

For the sake of illustration, here we focus on the role of internal noise stemming from the reactions (1), and ignore 
spatial degrees of freedom. As detailed in a forthcoming publication [18], and following the reasoning presented 
in [16], in a proper continuum limit, the spatial dispersal of individuals is accounted by diffusive terms in the SPDE (3). 

As in the main text, the overall number of individuals is denoted N and s = (a, b, c) stands for the frequencies (or 
densities) of the species A, B, and C in the population (i.e. s\ — a, s 2 = b and s 3 = c). The master equation giving 
the time-evolution of the probability P(s, t) of finding the system in the state s at time t then reads 

d t P{s, t)=J2 { p ( s + Ss > *)W( S + Ss -» a) - P{s, t)W(s -» s + Ss)} , (8) 

Ss 

where W(s — > s + Ss) denotes the transition probability from state s to the state s + Ss within one time step (loss 
term), W(s + Ss — > s) is the analogous gain term, and the summation extends over all possible changes Ss. As an 
example, the relevant changes Ssi = Sa in the density a resulting from the basic reactions (1) are 8si = 5a = —1/N 
in the third reaction, Ss± = l/N in the fourth, and zero in all others. We also choose the unit of time such that, on 
average, every individual reacts once per time step. The transition rates resulting from the reactions (1) then read 
W = Naac for the reaction CA C(d (the prefactor of N enters due to our choice of time scale, where N reactions 
occur in one unit of time) and W = Nfxa{l — a — b — c) for A(d — ^ AA. Transition probabilities associated with all 
other reactions (1) follow similarly. 

The system size, or Kramers-Moyal, expansion (SZE) [16] of the Master equation is an expansion in the increment 
Ss, which is proportional to A -1 . Therefore, the SZE may be understood as an expansion in the inverse system size 
A -1 . To second order in Ss, it yields the (generic) Fokker- Planck equation [16]: 

d t P(s,t) = -ai[ai(a)P(8,t)] + ^d l d J [B lJ (s)P(s,t)} . (9) 

For the system under consideration, in the above the indices i,j G {1,2,3} and the summation convention in ([9|) 
implies sums carried over them. According to the Kramers-Moyal expansion (or SZE), the quantities ai and Bij are 
given by [16] 

cti (s) = 2^ 5siW(s — > s + Ss) , 

Ss 

B t3 {s) =^2s Sl S Sj W(s -^s + Ss). (10) 

Ss 

Note that B is symmetric. For the sake of clarity, we outline the calculation of cti(s): The relevant changes Ss\ = Sa 
result from the third and fourth reactions in (1), as described above. The corresponding rates respectively read 
W = Naac and W = Nfia(l — a — b — c), resulting in ai(s) = fia(l — a — b — c) — aac. The other quantities are 
computed analogously. All explicit expressions for a»(s) and Bij(s) will be derived and given in detail in [18]. The 
well-known correspondence between Fokkcr-Planck equations and Ito calculus [16] implies that © is equivalent to 
the following set of Ito stochastic differential equations (with the above summation convention) : 

d t a = oti +Cij£j , 
d t b = a 2 +C 2 j£j , 

d t c = a 3 + C 3j £ j , (11) 

where the matrix C is defined from B via the relation CC T = B [16], and the &'s denote (uncorrelated) Gaussian 
white noise terms. Note that for the model under consideration B is diagonal and one can therefore always choose a 
diagonal matrix for C [see Eq. (5)], with only C^s contributing to the right-hand side of Eqs. (11). 
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In Rcf. [18], wc demonstrate that (in a proper continuum limit) spatial degrees of freedom and exchange processes 
simply yield additional diffusive terms V 2 s^ in (fTTj) . This leads to the SPDE (3), given and discussed in the main 
text, in which the £ 2 's still denote Gaussian white noise contributions. 

GENERALIZATION IN ARBITRARY DIMENSIONS 

The SPDE (3) and the CGLE (6), as well as the analytical predictions (7), are valid in all spatial dimensions. While 
in the text, for the sake of specificity, we mainly focus on the two-dimensional situation, here we comment on some 
features of the one-dimensional version of the system, as well as on some properties of the general situation in higher 
dimensions. 

In one dimension and in the absence of exchange processes (mixing), our model is expected to exhibit coarsening 
like the cyclic Lotka-Volterra model (see e.g. [12]). The same phenomenon still occurs in the presence of very slow 
mixing rate (for the model under consideration, this means very low values of the exchange rate e). On the other 
hand, and as shown in the main text, in the presence of (finite) mixing the systems' behavior is aptly described by 
the SPDE (3). The underlying CGLE (6) predicts the propagation of traveling waves, with velocity and wavelength 
still given by the analytical expression (7). Similarly to what has been found in the two-dimensional system [6], if 
the exchange rate is "moderate", i.e. below a certain mobility threshold but still finite [6,15], the system confirms 
these predictions and the species coexist. However, due to stochastic events and the presence of absorbing boundaries, 
some domains will occasionally merge, resulting in growing domains. This coarsening phenomenon will happen on a 
much longer time-scale than in the absence of the mixing processes. Above the threshold value for the diffusivity, the 
system is well-mixed, the underlying spatial structure plays no role, and the description in terms of the rate equations 
(2) is valid. In this mean-field scenario, the system approaches an absorbing steady state and no patterns emerge. 

In higher dimensions (i.e. in dimensions d > 2), the description of the stochastic spatial system in terms of the 
SPDE (3) and of the CGLE (7) is again valid for "moderate" (or intermediate) mixing (i.e. for a finite mobility rate 
which is below a certain critical threshold, see [6]). In this situation, and as discussed in the main text, the description 
in terms of the CGLE (6) is qualitatively valid and predicts the emergence of moving spiral waves in two dimensions 
(which is the case discussed in detail in the main text, see also [6,15]), and to "scroll waves", i.e. vortex filaments, in 
three dimensions [24]. 



